rm(list=ls())
library(tidyverse)
library(ggpubr)
library(haven)

anes <- read_spss("ANES files/anes_timeseries_cdf_spss_20220916.sav")


anes2 <- anes %>%
  select(c(weight_full = VCF0009z, #weights for full sample
           #weight_mode = VCF0009x, # FTF weights DROP 2020
           #weight_web = VCF0009y, # web weights DROP 2020
           weight_post_full = VCF9999, # post-election full weights
           year = VCF0004, # year
           FT_dem = VCF0218 , # Dem Party FT
           FT_rep = VCF0224, # Rep Party FT
           pid = VCF0301, # PID
           mode = VCF0017, #mode of survey, only want 0 (FTF) or 4 (web) later
           extEfcacy = VCF0648,
           satDem = VCF9255
  )) %>%
  subset(pid != 4 & pid > 0)%>% #drop pure independents and NAs
  sapply(as.numeric) %>% data.frame()


A24 <-  read_spss("ANES files/anes_timeseries_2024_spss_20250430.sav")

anes24 <- A24 %>% #2024 ANES
  mutate(year = 2024) %>%     
  select(weight_post_full =  V240107b,   # weights, Post full ANES sample (FTF + panel + web + PAPI)
         weight_full =   V240107a, # 
         #weight_ftf = V240101a, # FTF weights DROP 2020
         #weight_web = V240102a, # web weights DROP 2020
         FT_dem = V241166,
         FT_rep = V241167,
         pid =  V241227x,
         mode = V240002a, #only want 1 (FTF) or 2 (web)
         year,
         satDem = V242439,
         extEfcacy_official_dont_care = V242200, # POST: [STD] PUBLIC OFFICIALS DON’T CARE WHAT PEOPLE THINK, USE IT TO BUILD EFFICACY INDEX. 5-PT SCALE HERE UNLIKE ANES CDF
         extEfcacy_have_no_say = V242201 # POST: [STD] HAVE NO SAY ABOUT WHAT GOVERNMENT DOES, USE IT TO BUILD EFFICACY INDEX. 5-PT SCALE HERE UNLIKE ANES CDF     
  )%>%
  subset(pid != 4 & pid > 0)%>% #drop pure independents and NAs
  subset(FT_dem >= 0 & FT_rep >=0)%>%
  sapply(as.numeric) %>% data.frame() %>%
  mutate(satDem = as.numeric(recode(satDem,
                                    "1" = "1",
                                    "2" = '2',
                                    "4" = "3",
                                    "5" = "4"
  )))
 # mutate(extEfcacy = (extEfcacy_official_dont_care+extEfcacy_have_no_say)/2)

anes3 <- bind_rows(anes2, anes24)

#check
summary(anes3$year) #perf

##drop FT NAs

##########################################################################
#create pid dummy
anes3$pid_dum <- ifelse(anes3$pid < 4, "Democrat", "Republican")

#create inparty and outparty
anes3$inparty <- ifelse(anes3$pid_dum=="Democrat", "Democrat", "Republican")
anes3$outparty <- ifelse(anes3$pid_dum=="Democrat", "Republican", "Democrat")

#link FTs to in and out_party
anes3$inparty_FT <- ifelse(anes3$pid_dum=="Democrat", anes3$FT_dem, anes3$FT_rep)
anes3$outparty_FT <- ifelse(anes3$pid_dum=="Democrat", anes3$FT_rep, anes3$FT_dem)


###caclulate AP as IN_FT - OUT_FT
anes3$AP <- anes3$inparty_FT - anes3$outparty_FT


####
polar <- anes3 %>%
  select(year, AP, pid_dum, weight_full) %>%
  drop_na() %>%
  group_by(year) %>%
  summarise(value = weighted.mean(AP, weight_full),
            sd = sd(AP*weight_full),
            se = sd/sqrt(n()))


b<- ggplot(polar, aes(x=year, y=value)) +
  geom_line(linetype="dashed") +
  geom_point() +
  #  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin = (value -1.96*se), 
                    ymax = (value + 1.96*se)), width = 0.4) +
  theme_bw() +
  labs(x="", y="Mean",
       title="B",
       subtitle="Feeling thermometer ratings of one's own party (0 to 100) minus 
Feeling thermometer ratings of one’s rival party (0 to 100). Only partisans 
(including partisan leaners). N = 43,062 (ANES).
Means and 95% Confidence Intervals.")+
  scale_y_continuous(limits=c(0,100))

b

###satdem
satis <- anes3 %>%
  select(year, weight_full, satDem) %>%
  drop_na() %>%
  subset(satDem >=0) %>%
  mutate(satDem_reversed = ((satDem - 4)*-1)+1) %>%
  group_by(year)%>%
  summarise(value = weighted.mean(satDem_reversed, weight_full),
            sd = sd(satDem_reversed*weight_full),
            se = sd/sqrt(n()))
####CHECK


a<- ggplot(satis, aes(x=year, y=value)) +
  geom_point() +
  geom_line(linetype="dashed")+
  geom_errorbar(aes(ymin = (value -1.96*se), 
                    ymax = (value + 1.96*se)), width = 0.4) +
  theme_bw() +
  scale_x_continuous(breaks=c(2005,2010,2015,2020, 2024))+
  labs(x="", y="Category",
       title="A",
       subtitle="On the whole, are you satisfied with the way
democracy works in the United States? N=21,406 (ANES)
Means and 95% Confidence Intervals") +
  scale_y_continuous(limits=c(1,4),labels = c("Not at all satisfied",
                                              "Not very satisfied",
                                              "Fairly satisfied",
                                              "Very Satisfied"))
a

eff <- anes3 %>%
  subset((extEfcacy_official_dont_care >=0 & extEfcacy_have_no_say>=0) |
           (is.na(extEfcacy_official_dont_care) & is.na(extEfcacy_have_no_say)))


eff$extEfcacy <- ifelse(eff$year==2024, ((eff$extEfcacy_official_dont_care + eff$extEfcacy_have_no_say)/2)*10, eff$extEfcacy)

eff <- eff %>%
  select(year, weight_full, extEfcacy) %>%
  drop_na() %>%
  group_by(year) %>%
  summarise(value = weighted.mean(extEfcacy, weight_full),
            sd = sd(extEfcacy*weight_full),
            se = sd/sqrt(n()))

c <- ggplot(eff, aes(x=year, y=value)) +
  geom_line(linetype="dashed") +
  geom_point() +
  #  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin = (value -1.96*se), 
                    ymax = (value + 1.96*se)), width = 0.4) +
  theme_bw() +
  scale_y_continuous(limits=c(0,100)) +
  labs(x="", y="Mean",
       title="C",
       subtitle="Based on two items: People like me don't have any say
about what the government does, Agree (0), Neither (50), Disagree (100) /
I don't think public officials care much what people like me think, 
Agree (0), Neither (50), Disagree (100). N = 59,492 (ANES). 
Means and 95% Confidence Intervals.")

png(file = "PNAS/panel2.png",
    height = 6,
    width = 18,
    units = "in",       # Use inches to match pdf dimensions
    res = 300)          # Set resolution (DPI) for high-quality output

print(ggarrange(a,b,c,
                nrow = 1, ncol = 3, common.legend = T, legend = "bottom"))  

dev.off()

################### analyses
satdem2008 <- anes3 %>%
  subset(year==2008) %>%
  mutate(satDem_reversed = ((satDem - 4)*-1)+1) %>%
  select(satDem_reversed, weight_full) %>%
  summarise(satDem_reversed=satDem_reversed*weight_full) %>%
  drop_na()

satdem2012 <- anes3 %>%
  subset(year==2012) %>%
  mutate(satDem_reversed = ((satDem - 4)*-1)+1) %>%
  select(satDem_reversed, weight_full) %>%
  summarise(satDem_reversed=satDem_reversed*weight_full) %>%
  drop_na()

t.test(satdem2008, satdem2012,
       conf.level = 0.95)


sd(satdem2008$satDem_reversed)

sd(satdem2012$satDem_reversed)


eff <- anes3 %>%
  subset((extEfcacy_official_dont_care >=0 & extEfcacy_have_no_say>=0) |
           (is.na(extEfcacy_official_dont_care) & is.na(extEfcacy_have_no_say)))


eff$extEfcacy <- ifelse(eff$year==2024, ((eff$extEfcacy_official_dont_care + eff$extEfcacy_have_no_say)/2)*10, eff$extEfcacy)



eff00 <- eff %>%
  subset(year==2000) %>%
  mutate(extEfcacy = extEfcacy*weight_full)%>%
  select(extEfcacy) %>%
  drop_na()

eff24 <- eff %>%
  subset(year==2024) %>%
  mutate(extEfcacy = extEfcacy*weight_full)%>%
  select(extEfcacy) %>%
  drop_na()

sd(eff24$extEfcacy)
sd(eff00$extEfcacy)

mean(eff24$extEfcacy)
mean(eff00$extEfcacy)

t.test(eff00, eff24,
       conf.level = 0.95)


eff52 <- eff %>%
  subset(year==1952) %>%
  mutate(extEfcacy = extEfcacy*weight_full)%>%
  select(extEfcacy) %>%
  drop_na()

sd(eff52$extEfcacy)
mean(eff52$extEfcacy)



outpart00 <- anes3 %>%
  subset(year==2000) %>%
  select(outparty_score = outparty_FT, weight_full) %>%
  mutate(outparty_score = outparty_score*weight_full)%>%
  drop_na()

outpart24 <- anes3 %>%
  subset(year==2024) %>%
  select(outparty_score = outparty_FT, weight_full) %>%
  mutate(outparty_score = outparty_score*weight_full)%>%
  drop_na()

sd(outpart00$outparty_score)
sd(outpart24$outparty_score)


t.test(outpart00$outparty_score, outpart24$outparty_score,
       conf.level = 0.95)





### Appendix Mode Plot
modeplot <- anes3 %>%
  subset((year==2024 & mode %in% c(1,2)) | (year!= 2024 & mode %in% c(0,4)))%>%
  select(year, weight_full, mode, AP) %>%
  mutate(FTF = as.character(recode(mode,
                                   `0` ="FTF",
                                   `1` ="FTF",
                                   `2` ="Web",
                                   `4` = "Web")))


mode_summary <- modeplot%>%
  filter(!is.na(AP)) %>%
  group_by(year, FTF) %>%
  summarise(
    weighted_mean = weighted.mean(AP, weight_full),
    sd = sd(AP*weight_full),
    se= sd/sqrt(n()),
    .groups = "drop"
  )

ggplot(mode_summary, aes(x = year, y = weighted_mean, color = FTF, group = FTF)) +
  geom_line() +
  geom_point() +
    theme_bw() +
  theme(legend.position = "bottom")+
  geom_errorbar(aes(ymin = weighted_mean-se*1.96, ymax = weighted_mean+se*1.96), width=0.5) +
  labs(title = "Polarization across Survey Modes",
       x = "Year",
       y = "Mean",
       color = "Survey Mode") +
    scale_y_continuous(limits=c(0,100))
ggsave("PNAS/appendixB1.png")


